/* M. Nielsen Sept 2008 mniel@cbs.dtu.dk */

/* 
Copyright (C) 2008-2015 Danish Technical University

This suite of programs and library routine is free software. You can 
redistribute it and/or modify it under the terms of the GNU General Public 
License as published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

In other words, you are free to modify, copy, or redistribute this
source code and its documentation in any way you like, but you must
distribute all derivative versions as free software under the same
terms that I've provided my code to you (i.e. the GNU General Public
License). This precludes any use of the code in proprietary or
commercial software unless your source code is made freely available.

If you wish to use the code under a different Open Source license
that's not compatible with the GPL (like the Artistic License, BSD
license, or the Netscape Public License), please contact me
(Morten Nielsen, mniel@cbs.dtu.dk) for permission.

Incorporation into commercial software under non-GPL terms is possible
by obtaining a specially licensed version from The Danish Technical University.
Contact Morten Nielsen (mniel@cbs.dtu.dk) to arrange licensing terms.

This software is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
*/


#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "utils.h"

FILENAME        p_blrealmat;
int	p_verbose;
float	p_beta;
int	p_seqw;

PARAM   param[] = {
        "-blf", VFNAME p_blrealmat, "Real value blosum matrix filename", "$HOME/UBA/algobio/algobio/data/blosum62.freq_rownorm",
	"-v", VSWITCH	p_verbose, "Verbose mode", "0",
	"-b", VFLOAT	p_beta, "Weight on prior", "50.0",
	"-sw", VSWITCH	p_seqw, "Include sequece weighting", "0",
	0
};

float   *init_bg( )

{
        float   *bg;
       
        bg = fvector(0, 19 );
       
        bg[0]  = 0.074;
        bg[1]  = 0.052;
        bg[2]  = 0.045;
        bg[3]  = 0.054;
        bg[4]  = 0.025;
        bg[5]  = 0.034;
        bg[6]  = 0.054;
        bg[7]  = 0.074;
        bg[8]  = 0.026;
        bg[9]  = 0.068;
        bg[10] = 0.099;
        bg[11] = 0.058;
        bg[12] = 0.025;
        bg[13] = 0.047;
        bg[14] = 0.039;
        bg[15] = 0.057;
        bg[16] = 0.051;
        bg[17] = 0.013;
        bg[18] = 0.032;
        bg[19] = 0.073;

        return( bg );
}

float	sequence_weighting( PEPLIST *peplist, char *alphabet, int alen )

{
	PEPLIST	*pl;
	int	*r, **s;
	int	len;
	float	nseq;
	int	i,j, iv;
	float	w;

	len = peplist->len;

	r = ivector( 0, len-1 );
	s = imatrix( 0, len-1, 0, alen-1 );

	/* s is a matrix containing the number occurrence of each amino acids at each position */
	for ( pl=peplist; pl; pl=pl->next ) {

		for ( i=0; i<len; i++ ) {

			iv = pl->iv[i];

			s[i][iv]++;

		}

	}

	/* r is a vector containing the number of different amino acids at each position */
	for ( i=0; i<len; i++ ) {

		for ( j=0; j<alen; j++ ) {

			if( s[i][j] > 0)
				r[i]++;

		}
	}

	for ( pl=peplist; pl; pl=pl->next ) {

		pl->score = 0.0;

		for ( i=0; i<len; i++ ) {

			iv = pl->iv[i];

			w = 1.0/(r[i]*s[i][iv]);
			pl->score += w;
		}

	}

	nseq = 0;
	for ( i=0; i<len; i++ ) 
		nseq += r[i];

	nseq /= len;

	return( nseq );
}

main( int argc, char *argv[] )

{
	PEPLIST	*peplist, *pl;
	float		*bg;
	char		*bl_alphabet;
	float		nseq;
	float		**aacounts;
	float		**aafrec;
	float		**aalo;
	float		**blm;
	int		i, j, k, ix;
	int		len, alen, c;
	float		g, sum;

	pparse( &argc, &argv, param, 1, "peptide_file" );

	/* The peplist_read routine reads in a linked list of PEPLIST elements from a given file */
	/* The PEPLIST elements are defined in utils and contain the peptide sequence, length, score, ... */
	peplist = peplist_read( argv[1] );
	if ( peplist == NULL ) {
			printf( "No Peptides sequences read from file %s\n", argv[1] );
			exit( 0 );
	}

	/* Check that all peptides have equal length */
	for ( nseq=0, pl=peplist; pl; pl=pl->next, nseq+=1.0 ) {

		if ( pl->len != peplist->len ) {

			printf( "Inconsistent peptide length %i %s %i %s. Exit\n", 
				pl->len, pl->pep, peplist->len, peplist->pep );
			exit( 1 );
		}
	}

	/* Read real value Blosum substitution frequence matrix */
	/* The read_realblosum subroutine reads a blosum substitution frequency matrix, and the corresponding
	alphabet (the amino acids order) from a given file. The matrix is stored in the variable blm and the alphabet
	in the variable alphabet (a string) */

	blm = read_realblosum( p_blrealmat, &bl_alphabet );

	/* Assign peptide amino acid position in the Blosum alphabet to the iv vector for each peptide element */
	for ( pl=peplist; pl; pl=pl->next ) {

		/* Allocate an integer vector (iv) to store the translation from amino acids to integers */
		pl->iv = ivector( 0, pl->len-1 );
		pl->score = 1.0;

		for ( i=0; i<pl->len; i++ ) {
			c = strpos( bl_alphabet, (pl->pep)[i] );

			if ( c < 0 ) {
				printf( "Error. Peptide AA %c not in alphabet %s. Exit\n", (pl->pep)[i], bl_alphabet );
				exit( 1 );
			}

			pl->iv[i] = c;
		}
	}

	/* Read Background frequencies */

	bg = init_bg( );

	len = peplist->len;
	alen = strlen( bl_alphabet );

	/* Allocate matrices for calculation */

	aacounts = fmatrix(0, len-1, 0, alen-1);
	aafrec = fmatrix(0, len-1, 0, alen-1);
	aalo = fmatrix(0, len-1, 0, alen-1);

	/* if option -sw is switched on call sequence_weighting subroutine */
	/* This routine assigns the weight of each peptide to the peplist variable score ie. pl->score */
	/* It returns the value for alpha to be used when calculating the combined amino acids frequencies */

	if ( p_seqw )
		nseq = sequence_weighting( peplist, bl_alphabet, alen );

	printf( "# Number of sequences nseq %f\n", nseq );

	/* Now find the amino acid counts at each peptide position */
	for ( pl=peplist; pl; pl=pl->next ) {

		for ( i=0; i<len; i++ ) {

			ix = pl->iv[i];

			aacounts[i][ix] += pl->score;

		}

	}

	/* Find sum of sequence weights */
	sum = 0.0;
	for ( pl=peplist; pl; pl=pl->next )
		sum += pl->score;

	/* Normalize the amino acids counts to frequencies */
	for ( i=0; i<len; i++ ) {

		for ( j=0; j<alen; j++ )
			aacounts[i][j] /= ( sum > 0 ? sum : 1.0 );
	}

	/* Calculate pseudo frequency g, and combined frequence (a*f+b*g)/(a+b), where a is nseq-1 and b is 
	(the command line parameter) weight on prior */
	float a = nseq - 1;
	for ( i=0; i<len; i++ ) {

		for ( j=0; j<alen; j++ ) {

			/* g_j = sum_k f_k * q(j|k) */
			/* f_k =  aacounts[i][k]; */
			/* q(j|k) = blm[k][j] */

			g = 0.0;
			
			for( k=0; k<alen; k++)
				g += aacounts[i][k] * blm[k][j];
			float_t	 f = aacounts[i][j];
			aafrec[i][j] = (a*f+p_beta*g)/(a+p_beta);

		}

	}

	/* Make log-odds scores */
	/* aalo[i][j] = 2 * log ( aafreq/bg ) / log( 2.0 ) */
		
	for ( i=0; i<len; i++ ) 
		for ( j=0; j<alen; j++ ) 
			if ( aafrec[i][j] > 0 )
				aalo[i][j] = 2*log( aafrec[i][j] / bg[j] )/log(2.0);

	/* 
	*	Print logodds  matrix
	*/

	printf( "\nLast position-specific scoring matrix computed       \n");

	printf( "      " );
	for (i=0;i<alen;i++)
		printf( "%7c ", bl_alphabet[i] );
	printf( "\n");

	for (i=0;i<len;i++){
		printf( "%3d A ",i+1);
		for (j=0;j<alen;j++)
			printf( "%7.3f ", aalo[i][j]);
		printf( "\n");
	}

	exit( 0 );
}
